rm(list = ls())
library(readstata13)
library(plyr)
library(plotrix)
library(foreign)
library(stargazer)
library(ggplot2)
library(stm)
library(dplyr)
library(zoo)
library(car)
library(plyr)
library(knitr)
library(qwraps2)
library(knitr)
library(xtable)
library(qwraps2)
library(tibble)
library(gdata)
library(tidyr)
library(data.table)
library(stringr)
library(reshape2)
library(ggrepel)
library(dfadjust)
library(lmtest)
library(summarytools)
library(kableExtra)
library(experiment)
library(AER)
library(sandwich)
#library(ivpack)
library(psych)
library(remotes)
library(stats)
library(ggpubr)

setwd("~/Dropbox/Participacion Politica FARC/Data_workshops/Replication/Dataverse files")

load("data_notimputed.RData")
load("data_imputed.RData")

imputedT <- subset(imputed, treatment=="1")
imputedT <- as_tibble(imputedT)
imputedT <- merge(imputedT, imputedT %>% dplyr::count(schooling), by = "schooling")
colnames(imputedT)[length(colnames(imputedT))] <- "N"

imputedC <- subset(imputed, treatment=="0")
imputedC <- as_tibble(imputedC)
imputedC <- merge(imputedC, imputedC %>% dplyr::count(schooling), by = "schooling")
colnames(imputedC)[length(colnames(imputedC))] <- "N"

##### Figure 3: Education and moderation before and after the workshop ##### 
c <- imputedC %>% group_by(schooling) %>% summarize(mean_trust_byestudio = mean(trust_score, na.rm = T))
c <- merge(c, imputedC %>% count(schooling), by = "schooling")
colnames(c)[3] <- "N"
c <- merge(c,imputedC %>% group_by(schooling) %>% summarize(mean_democracy_byestudio = mean(democracy_score, na.rm = T)))
c <- merge(c,imputedC %>% group_by(schooling) %>% summarize(mean_moderation_byestudio = mean(moderation_score, na.rm = T)))
c <- merge(c,imputedC %>% group_by(schooling) %>% summarize(mean_ideomoderation_byestudio = mean(moderation_platform, na.rm = T)))
c <- merge(c,imputedC %>% group_by(schooling) %>% summarize(mean_behmoderation_byestudio = mean(moderation_alliance, na.rm = T)))
c <- merge(c,imputedC %>% group_by(schooling) %>% summarize(mean_ideology_byestudio = mean(ideology_personal, na.rm = T)))
t <- imputedT %>% group_by(schooling) %>% summarize(mean_trust_byestudio = mean(trust_score, na.rm = T))
t <- merge(t, imputedT %>% count(schooling), by = "schooling")
colnames(t)[3] <- "N"
t <- merge(t,imputedT %>% group_by(schooling) %>% summarize(mean_democracy_byestudio = mean(democracy_score, na.rm = T)))
t <- merge(t,imputedT %>% group_by(schooling) %>% summarize(mean_moderation_byestudio = mean(moderation_score, na.rm = T)))
t <- merge(t,imputedT %>% group_by(schooling) %>% summarize(mean_ideomoderation_byestudio = mean(moderation_platform, na.rm = T)))
t <- merge(t,imputedT %>% group_by(schooling) %>% summarize(mean_behmoderation_byestudio = mean(moderation_alliance, na.rm = T)))
t <- merge(t,imputedT %>% group_by(schooling) %>% summarize(mean_ideology_byestudio = mean(ideology_personal, na.rm = T)))

f3a <- ggplot(data=c,
       aes(x=schooling, y=mean_moderation_byestudio)) +
  geom_point(aes(size = N), colour = "red")+ 
  geom_smooth(aes(weight = N), method='lm', se = F, size = 0.5, color = "black", na.rm = T) +
  labs(#title = "Ex-combatant education and moderation score",
       subtitle = "Control group (N = 117)",
       x = "Level of education",
       y = "Moderation score")+ 
  scale_y_continuous(limits = c(-2,1))+ 
  theme(text = element_text(size = 8))+
  theme(panel.background = element_blank())+
  theme(legend.position="none")+
  theme(axis.line = element_line(colour = "black")) 

f3b <- ggplot(data=t,
       aes(x=schooling, y=mean_moderation_byestudio)) +
  geom_point(aes(size = N), colour = "blue")+ 
  geom_smooth(aes(weight = N), method='lm', se = F, size = 0.5, color = "black", na.rm = T) +
  labs(#title = "Ex-combatant education and moderation score",
       subtitle = "Treatment group (N = 158)",
       x = "Level of education",
       y = "Moderation score")+ 
  scale_y_continuous(limits = c(-2,1))+ 
  theme(text = element_text(size = 8))+
  theme(panel.background = element_blank())+
  theme(legend.position="none")+
 theme(axis.line = element_line(colour = "black")) 

f3 <- ggarrange(f3a,f3b, legend = "none")

ggsave("f3.tiff",plot = f3, dpi = 300, width = 4, height = 1.75, unit = "in")
ggsave("f3.png",plot = f3, width = 4, height = 1.75, unit = "in")

##### Figure 4: Mediation of ideology on moderation score #####
library(multilevel)
library(mediation)
set.seed(2020)
imputed$treated <- as.factor(imputed$treatment)
subfarc <- imputed %>% dplyr::select(ideology_personal, treated, treatment,moderation_score)
subfarc$treated <- as.factor(subfarc$treatment)
subfarc <- na.omit(subfarc)
fitM <- lm(ideology_personal ~ treated, data = subfarc)
fitY <- lm(moderation_score ~ treated + ideology_personal, data = subfarc)
fitMed <- mediation::mediate(fitM, fitY, treat="treated", mediator="ideology_personal",boot = TRUE, sims = 500)
summary(fitMed)
png(file="f4.png")
plot(fitMed, main="Mediation of ideology of effect on moderation effect")
dev.off()

##### Figure 5: Education and moderation before and after the workshop ##### 
c <- imputedC %>% group_by(conflict_score) %>% summarize(mean_trust_byestudio = mean(trust_score, na.rm = T))
c <- merge(c, imputedC %>% count(conflict_score), by = "conflict_score")
colnames(c)[3] <- "N"
c <- merge(c,imputedC %>% group_by(conflict_score) %>% summarize(mean_democracy_byestudio = mean(democracy_score, na.rm = T)))
c <- merge(c,imputedC %>% group_by(conflict_score) %>% summarize(mean_moderation_byestudio = mean(moderation_score, na.rm = T)))
t <- imputedT %>% group_by(conflict_score) %>% summarize(mean_trust_byestudio = mean(trust_score, na.rm = T))
t <- merge(t, imputedT %>% count(conflict_score), by = "conflict_score")
colnames(t)[3] <- "N"
t <- merge(t,imputedT %>% group_by(conflict_score) %>% summarize(mean_democracy_byestudio = mean(democracy_score, na.rm = T)))
t <- merge(t,imputedT %>% group_by(conflict_score) %>% summarize(mean_moderation_byestudio = mean(moderation_score, na.rm = T)))


f5a <- ggplot(data=c,
       aes(x=conflict_score, y=mean_moderation_byestudio)) +
  geom_point(aes(size = N), colour = "red")+ 
  geom_smooth(aes(weight = N), method='lm', se = F, size = 0.5, color = "black", na.rm = T) +
  labs(#title = "Ex-combatant conflict experience and moderation score",
       subtitle = "Control group (N = 117)",
       x = "Conflict experience score",
       y = "Moderation score")+ 
  scale_y_continuous(limits = c(-2, 1))+ 
  scale_x_continuous(limits = c(-2, 1))+
  theme(axis.line = element_line(colour = "black")) +
  theme(text = element_text(size = 8))+
  theme(panel.background = element_blank())+
  theme(legend.position="none")

f5b <- ggplot(data=t,
       aes(x=conflict_score, y=mean_moderation_byestudio)) +
  geom_point(aes(size = N), colour = "blue")+ 
  geom_smooth(aes(weight = N), method='lm', se = F, size = 0.5, color = "black", na.rm = T) +
  labs(#title = "Ex-combatant conflict experience and moderation score",
       subtitle = "Treatment group (N = 158)",
       x = "Conflict experience score",
       y = "Moderation score")+ 
  scale_y_continuous(limits = c(-2, 1))+ 
  scale_x_continuous(limits = c(-2, 1))+
  theme(axis.line = element_line(colour = "black")) +
  theme(text = element_text(size = 8))+
  theme(panel.background = element_blank())+
  theme(legend.position="none")

f5 <- ggarrange(f5a,f5b, legend = "none")

ggsave("f5.tiff",plot = f5, dpi = 300, width = 4, height = 1.75, unit = "in")
ggsave("f5.png",plot = f5, width = 4, height = 1.75, unit = "in")

#### Figure A 1: Balance across demographic characteristics ####
library(ggplot2)
library(plyr)

age <- ggplot(farc, aes(x=age, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5) + 
  labs(title = " ",
       x = "Age",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA1_balance_age.png",plot = age)


race <- ggplot(farc, aes(x=race, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5) + 
  labs(title = " ",
       x = "Race",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA1_balance_race.png",plot = race)

schooling <- ggplot(farc, aes(x=schooling, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5) + 
  labs(title = " ",
       x = "Education",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA1_balance_schooling.png",plot = schooling)

pdr <- ggplot(farc, aes(x=pdr, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "Currently reincorporation process (PDR)",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA1_balance_pdr.png",plot = pdr)


voted <- ggplot(farc, aes(x=voted2016, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "Voted in 2016",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA1_balance_voted.png",plot = voted)

#### Figure A 2: Balance across conflict experience  ####

entry_year <-  ggplot(farc, aes(x=entry_year, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "Year of entry into FARC (guerrilla)",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_yearentry.png",plot = entry_year)

entry_age <- ggplot(farc, aes(x=entry_age, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5) + 
  labs(title = " ",
       x = "Age of entry into FARC (guerrilla)",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_ageentry.png",plot = entry_age)


yearsin_farc <- ggplot(farc, aes(x=yearsin_farc, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "Years spent with FARC (guerrilla)",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_timespent.png",plot = yearsin_farc)


ideological_education <- ggplot(farc, aes(x=ideological_education, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "During conflict, received ideological training",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_ideological_education.png",plot =ideological_education)


commander_interest <- ggplot(farc, aes(x=commander_interest, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "During conflict, commander cared about what respondent had to say",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_commander_interest.png",plot =commander_interest)


debate <- ggplot(farc, aes(x=debate, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "During conflict, debated politics with unit",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_debate.png",plot =debate)


shot <- ggplot(farc, aes(x=shot, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "During conflict, was shot at",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_shot.png",plot =shot)


lifethreat <- ggplot(farc, aes(x=lifethreat, fill = treated, color=treated)) +
  geom_histogram(position="identity", alpha = 0.5)+ 
  labs(title = " ",
       x = "During conflict, life was threatened",
       y = "Frequency")+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA2_balance_lifethreat.png", plot =lifethreat)   
#### Figure A 5: Difference in means across outcomes (a-d) #### 
#preparing data for the plots
farcT <- subset(farc, treatment=="1")
farcC <- subset(farc, treatment=="0")

#### Figure A 5a Trust in democracy #### 
d1T <- mean(farcT$system_inclusive, na.rm = T)
d1C <- mean(farcC$system_inclusive, na.rm = T)
dse1T <- std.error(farcT$system_inclusive, na.rm = T)
dse1C <- std.error(farcC$system_inclusive, na.rm = T)

d2T <- mean(farcT$democracy_best, na.rm = T)
d2C <- mean(farcC$democracy_best, na.rm = T)
dse2T <- std.error(farcT$democracy_best, na.rm = T)
dse2C <- std.error(farcC$democracy_best, na.rm = T)

d3T <- mean(farcT$farc_goals, na.rm = T)
d3C <- mean(farcC$farc_goals, na.rm = T)
dse3T <- std.error(farcT$farc_goals, na.rm = T)
dse3C <- std.error(farcC$farc_goals, na.rm = T)

d4T <- mean(farcT$voice_ingov, na.rm = T)
d4C <- mean(farcC$voice_ingov, na.rm = T)
dse4T <- std.error(farcT$voice_ingov, na.rm = T)
dse4C <- std.error(farcC$voice_ingov, na.rm = T)

d5T <- mean(farcT$mecdem_efficient, na.rm = T)
d5C <- mean(farcC$mecdem_efficient, na.rm = T)
dse5T <- std.error(farcT$mecdem_efficient, na.rm = T)
dse5C <- std.error(farcC$mecdem_efficient, na.rm = T)

means3 <- c(d1C, d1T, d2C, d2T, d3C, d3T, d4C, d4T, d5C, d5T)
se3 <- c(dse1C, dse1T, dse2C, dse2T, dse3C, dse3T, dse4C, dse4T, dse5C, dse5T)

dataplotDEMO <- data.frame(x3 = c('Inclusive', 'Inclusive', 'Best', 'Best', 'FARC goals', 'FARC goals', 'Voice', 'Voice', 'Efficient', 'Efficient'),
                           means3, se3,
                           Treatment = c('Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment'))

democracy <- ggplot(dataplotDEMO, aes(x=x3, y=means3, fill=Treatment)) + 
  geom_bar(stat="identity",
           position=position_dodge()) +
  geom_errorbar(aes(ymin=means3-se3, ymax=means3+se3), width=.2,
                position=position_dodge(.9)) +
  theme(text = element_text(size = 15)) +
  scale_x_discrete(name ="Democracy")+
  scale_y_continuous(name = "Mean agreement", limits = c(0, 10))+
  theme(panel.background = element_blank())+
  theme_bw()
ggsave("fA5_a.png",plot = democracy)
#### Figure A 5b Trust in institutions #### 
#trust in government 
c1T <- mean(farcT$trust_gov, na.rm = T)
c1C <- mean(farcC$trust_gov, na.rm = T)
se1T <- std.error(farcT$trust_gov, na.rm = T)
se1C <- std.error(farcC$trust_gov, na.rm = T)
#trust in mayors office
c2T <- mean(farcT$trust_mayors, na.rm = T)
c2C <- mean(farcC$trust_mayors, na.rm = T)
se2T <- std.error(farcT$trust_mayors, na.rm = T)
se2C <- std.error(farcC$trust_mayors, na.rm = T)
#trust in municipal council
c3T <- mean(farcT$trust_council, na.rm = T)
c3C <- mean(farcC$trust_council, na.rm = T)
se3T <- std.error(farcT$trust_council, na.rm = T)
se3C <- std.error(farcC$trust_council, na.rm = T)
#trust in justice system
c4T <- mean(farcT$trust_justice, na.rm = T)
c4C <- mean(farcC$trust_justice, na.rm = T)
se4T <- std.error(farcT$trust_justice, na.rm = T)
se4C <- std.error(farcC$trust_justice, na.rm = T)
#trust in JEP 
c5T <- mean(farcT$trust_jep, na.rm = T)
c5C <- mean(farcC$trust_jep, na.rm = T)
se5T <- std.error(farcT$trust_jep, na.rm = T)
se5C <- std.error(farcC$trust_jep, na.rm = T)

means <- c(c1C, c1T, c2C, c2T, c3C, c3T, c4C, c4T, c5C, c5T)
se <- c(se1C, se1T, se2C, se2T, se3C, se3T, se4C, se4T, se5C, se5T)

dataplot <- data.frame(x = c('Government', 'Government', 'Mayor', 'Mayor', 'Consejo', 'Consejo', 'Justice', 'Justice', 'Justice Peace', 'Justice Peace'),
                       means, se,
                       Treatment = c('Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment'))

institutions <- ggplot(dataplot, aes(x=x, y=means, fill=Treatment)) + 
  geom_bar(stat="identity",
           position=position_dodge()) +
  geom_errorbar(aes(ymin=means-se, ymax=means+se), width=.2,
                position=position_dodge(.9)) +
  theme(text = element_text(size = 15)) +
  scale_x_discrete(name ="Trust in institutions")+
  scale_y_continuous(name = "Mean agreement", limits = c(0, 10))+
  theme(panel.background = element_blank())+
  theme_bw()
ggsave("fA5_b.png",plot = institutions)

#### Figure A 5c Participation in democracy #### 
a1T <- mean(farcT$plan_vote2019, na.rm = T)
a1C <- mean(farcC$plan_vote2019, na.rm = T)
ase1T <- std.error(farcT$plan_vote2019, na.rm = T)
ase1C <- std.error(farcC$plan_vote2019, na.rm = T)

a2T <- mean(farcT$support_farc_candidate, na.rm = T)
a2C <- mean(farcC$support_farc_candidate, na.rm = T)
ase2T <- std.error(farcT$support_farc_candidate, na.rm = T)
ase2C <- std.error(farcC$support_farc_candidate, na.rm = T)

a3T <- mean(farcT$plan_campaigning, na.rm = T)
a3C <- mean(farcC$plan_campaigning, na.rm = T)
ase3T <- std.error(farcT$plan_campaigning, na.rm = T)
ase3C <- std.error(farcC$plan_campaigning, na.rm = T)

means2 <- c(a1C, a1T, a2C, a2T, a3C, a3T)
se2 <- c(ase1C, ase1T, ase2C, ase2T, ase3C, ase3T)

dataplota <- data.frame(x2 = c('Vote', 'Vote', 'Support Candidate', 'Support Candidate', 'Campaign', 'Campaign'),
                        means2, se2,
                        Treatment = c('Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment'))

participation <- ggplot(dataplota, aes(x=x2, y=means2, fill=Treatment)) + 
  geom_bar(stat="identity",
           position=position_dodge()) +
  geom_errorbar(aes(ymin=means2-se2, ymax=means2+se2), width=.2,
                position=position_dodge(.9)) +
  theme(text = element_text(size = 15)) +
  scale_x_discrete(name ="Future political participation") +
  scale_y_continuous(name = "Mean participation (N/Y)", limits = c(0, 1))+
  theme(panel.background = element_blank())+
  theme_bw()
ggsave("fA5_c.png",plot = participation)
#### Figure A 5d Support for party moderation #### 
e7T <- mean(farcT$moderation_platform, na.rm = T)
e7C <- mean(farcC$moderation_platform, na.rm = T)
ese7T <- std.error(farcT$moderation_platform, na.rm = T)
ese7C <- std.error(farcC$moderation_platform, na.rm = T)

e8T <- mean(farcT$moderation_alliance, na.rm = T)
e8C <- mean(farcC$moderation_alliance, na.rm = T)
ese8T <- std.error(farcT$moderation_alliance, na.rm = T)
ese8C <- std.error(farcC$moderation_alliance, na.rm = T)

means4 <- c(e7C, e7T, e8C, e8T)
se4 <- c(ese7C, ese7T, ese8C, ese8T)


dataplotPRAG <- data.frame(x4 = c( 'Platform', 'Platform', 'Alliances', 'Alliances'),
                           means4, se4,
                           Treatment = c('Control', 'Treatment', 'Control', 'Treatment'))

moderation <- ggplot(dataplotPRAG, aes(x=x4, y=means4, fill=Treatment)) + 
  geom_bar(stat="identity",
           position=position_dodge()) +
  geom_errorbar(aes(ymin=means4-se4, ymax=means4+se4), width=.2,
                position=position_dodge(.9)) +
  theme(text = element_text(size = 15)) +
  scale_x_discrete(name ="Support for party moderation")+
  scale_y_continuous(name = "Mean agreement", limits = c(0, 10))+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA5_d.png",plot = moderation)

#### Figure A 6: Difference in means across measures of ideology #### 
f4T <- mean(farcT$ideology_farc, na.rm = T)
f4C <- mean(farcC$ideology_farc, na.rm = T)
fse4T <- std.error(farcT$ideology_farc, na.rm = T)
fse4C <- std.error(farcC$ideology_farc, na.rm = T)

f5T <- mean(farcT$ideology_idealfarc, na.rm = T)
f5C <- mean(farcC$ideology_idealfarc, na.rm = T)
fse5T <- std.error(farcT$ideology_idealfarc, na.rm = T)
fse5C <- std.error(farcC$ideology_idealfarc, na.rm = T)

f6T <- mean(farcT$ideology_farcelections, na.rm = T)
f6C <- mean(farcC$ideology_farcelections, na.rm = T)
fse6T <- std.error(farcT$ideology_farcelections, na.rm = T)
fse6C <- std.error(farcC$ideology_farcelections, na.rm = T)

f7T <- mean(farcT$ideology_personal, na.rm = T)
f7C <- mean(farcC$ideology_personal, na.rm = T)
fse7T <- std.error(farcT$ideology_personal, na.rm = T)
fse7C <- std.error(farcC$ideology_personal, na.rm = T)

means5 <- c(f4C, f4T, f5C, f5T, f6C, f6T, f7C, f7T)
se5 <- c(fse4C, fse4T, fse5C, fse5T, fse6C, fse6T, fse7C, fse7T)


dataplotIDEO <- data.frame(x5 = c('FARC', 'FARC', 'Ideal FARC', 'Ideal FARC', 'Election FARC', 'Election FARC', 'Personal', 'Personal'),
                           means5, se5,
                           Treatment = c('Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment', 'Control', 'Treatment'))

ideology <- ggplot(dataplotIDEO, aes(x=x5, y=means5, fill=Treatment)) + 
  geom_bar(stat="identity",
           position=position_dodge()) +
  geom_errorbar(aes(ymin=means5-se5, ymax=means5+se5), width=.2,
                position=position_dodge(.9)) +
  theme(text = element_text(size = 15)) +
  scale_x_discrete(name ="Ideology")+
  scale_y_continuous(name = "Mean ideology (R - L)", limits = c(0, 10))+
  theme(panel.background = element_blank())+
  theme_bw()

ggsave("fA6.png",plot = ideology)

#### Figure A 7: Trust in institutions and in democracy as mediators of moderation #### 
library(multilevel)
library(mediation)
set.seed(2020)

imputed$treated <- as.factor(imputed$treatment)
subfarc <- imputed %>% dplyr::select(trust_score, treatment,moderation_score, treated)
subfarc <- na.omit(subfarc)
fitM <- lm(trust_score ~ treated, data = subfarc)
fitY <- lm(moderation_score ~ treated + trust_score, data = subfarc)
fitMeda <- mediation::mediate(fitM, fitY, treat="treated", mediator="trust_score",boot = TRUE, sims = 500)
#summary(fitMed)
png(file="fA7a.png")
plot(fitMeda, main="Trust in inst. as a mediator of effect on moderation effect")
dev.off()

subfarc <- imputed %>% dplyr::select(democracy_score, treatment,moderation_score, treated)
subfarc <- na.omit(subfarc)
fitM <- lm(democracy_score ~ treated, data = subfarc)
fitY <- lm(democracy_score ~ treated + democracy_score, data = subfarc)
fitMedb <- mediation::mediate(fitM, fitY, treat="treated", mediator="democracy_score",boot = TRUE, sims = 500)
#summary(fitMed)
png(file="fA7b.png")
plot(fitMedb, main="Trust in democracy as a mediator of moderation effect")
dev.off()

#### Figure A 8: Effect size by sites #### 
imputed <- merge(imputed, dplyr::select(farc, ID, front, position), by = "ID")

#Creating front level data to create mesaure of fractionalization
#creating variables that capture unit fracitonalization and workshop duration 
fronts <- imputed %>% dplyr::select(front,municipio)
fronts$front <- trimws(fronts$front)
fronts <- fronts %>% filter(fronts$front != "NA")

#wsnobs is the number of fronts that are within each municipality 
fronts <-  merge(fronts, fronts %>% group_by(municipio) %>% summarise(wsnobs = length(front)))

#wnobs is the number of individuals that are within each municipality 
fronts <- merge(fronts, imputed %>% group_by(municipio) %>% summarise(wnobs = length(municipio)))

#snobs is the number of individuals that are within each municipality 
fronts <- merge(fronts, fronts  %>% group_by(front, municipio) %>% summarise(snobs = length(front)))
fronts <- unique(fronts)

#s = for each unit within each workshop site, 1 - fractionalization measure
fronts$suw_2 <- (fronts$snobs / fronts$wnobs)^2

#merging data 
fronts <- merge(fronts, fronts %>% dplyr::group_by(municipio) %>% dplyr::summarise(uf = 1 - sum(suw_2)))

imputed <- merge(imputed, unique(dplyr::select(fronts, municipio, uf)))
f <- unique(dplyr::select(fronts,municipio,uf))
colnames(f)[1] <- "Location"


#Estimating effect sizes by site 
imputed$coyaima <- ifelse(imputed$municipio == "Coyaima", 1,0)
#N = 29
Coyaima <- imputed %>% dplyr::filter(imputed$municipio == "Coyaima")
c1 <- coeftest(lm(trust_score ~ assignment, weight = Coyaima$weight, data = Coyaima),vcov = vcovHC, type = "HC0")
c2 <- coeftest(lm(democracy_score ~ assignment, weight = weight, data = Coyaima),vcov = vcovHC, type = "HC0")
c3 <- coeftest(lm(moderation_platform ~ assignment, weight = weight, data = Coyaima),vcov = vcovHC, type = "HC0")

imputed$planadas <- ifelse(imputed$municipio == "Planadas", 1,0)
#N = 26
Planadas <- imputed %>% dplyr::filter(imputed$municipio == "Planadas")
p1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Planadas),vcov = vcovHC, type = "HC0")
p2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Planadas),vcov = vcovHC, type = "HC0")
p3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Planadas),vcov = vcovHC, type = "HC0")

imputed$iconozo <- ifelse(imputed$municipio == "Iconozo", 1,0)
#N = 14
Iconozo <- imputed %>% dplyr::filter(imputed$municipio == "Iconozo")
i1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Iconozo),vcov = vcovHC, type = "HC0")
i2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Iconozo),vcov = vcovHC, type = "HC0")
i3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Iconozo),vcov = vcovHC, type = "HC0")
i4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = Iconozo),vcov = vcovHC, type = "HC0")

imputed$florencia <- ifelse(imputed$municipio == "Florencia", 1,0)
#N = 13
Florencia <- imputed %>% dplyr::filter(imputed$municipio == "Florencia")
f1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Florencia),vcov = vcovHC, type = "HC0")
f2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Florencia),vcov = vcovHC, type = "HC0")
f3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Florencia),vcov = vcovHC, type = "HC0")
f4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = Florencia),vcov = vcovHC, type = "HC0")

imputed$montañita <- ifelse(imputed$municipio == "Montañita", 1,0)
#N = 32
Montañita <- imputed %>% dplyr::filter(imputed$municipio == "Montañita")
m1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Montañita),vcov = vcovHC, type = "HC0")
m2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Montañita),vcov = vcovHC, type = "HC0")
m3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Montañita),vcov = vcovHC, type = "HC0")
m4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = Montañita),vcov = vcovHC, type = "HC0")

imputed$bogota <- ifelse(imputed$municipio == "Bogotá", 1,0)
#N = 16
Bogota <- imputed %>% dplyr::filter(imputed$municipio == "Bogotá")
b1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Bogota),vcov = vcovHC, type = "HC0")
b2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Bogota),vcov = vcovHC, type = "HC0")
b3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Bogota),vcov = vcovHC, type = "HC0")
b4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = Bogota),vcov = vcovHC, type = "HC0")

imputed$villavicencio <- ifelse(imputed$municipio == "Villavicencio", 1,0)
#N = 22
Villavicencio <- imputed %>% dplyr::filter(imputed$municipio == "Villavicencio")
v1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Villavicencio),vcov = vcovHC, type = "HC0")
v2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Villavicencio),vcov = vcovHC, type = "HC0")
v3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Villavicencio),vcov = vcovHC, type = "HC0")
v4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = Villavicencio),vcov = vcovHC, type = "HC0")

imputed$vistahermosa <- ifelse(imputed$municipio == "Vista Hermosa", 1,0)
#N = 91
VistaHermosa <- imputed %>% dplyr::filter(imputed$municipio == "Vista Hermosa")
vh1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = VistaHermosa),vcov = vcovHC, type = "HC0")
vh2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = VistaHermosa),vcov = vcovHC, type = "HC0")
vh3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = VistaHermosa),vcov = vcovHC, type = "HC0")
vh4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = VistaHermosa),vcov = vcovHC, type = "HC0")

imputed$arauquita <- ifelse(imputed$municipio == "Arauquita", 1,0)
#N = 32
Arauquita <- imputed %>% dplyr::filter(imputed$municipio == "Arauquita")
a1 <- coeftest(lm(trust_score ~ assignment , weight = weight, data = Arauquita),vcov = vcovHC, type = "HC0")
a2 <- coeftest(lm(democracy_score ~ assignment , weight = weight, data = Arauquita),vcov = vcovHC, type = "HC0")
a3 <- coeftest(lm(moderation_platform ~ assignment , weight = weight, data = Arauquita),vcov = vcovHC, type = "HC0")
a4 <- coeftest(lm(moderation_alliance ~ assignment , weight = weight, data = Arauquita),vcov = vcovHC, type = "HC0")

stargazer(a1,b1,c1,f1,i1,m1,p1,v1, vh1)
location <- rbind.data.frame(as.data.frame(unique(imputed$municipio)),
                             as.data.frame(unique(imputed$municipio)))

colnames(location)[1] <- "Location"
location$Location <- sort(location$Location)
location$beta <- c("Effect","Constant","Effect","Constant","Effect","Constant","Effect","Constant","Effect","Constant","Effect","Constant","Effect","Constant","Effect","Constant","Effect","Constant")

location$trust <-  c(a1[2], a1[1], 
                     b1[2], b1[1],
                     c1[2], c1[1],
                     f1[2], f1[1],
                     i1[2], i1[1], 
                     m1[2], m1[1], 
                     p1[2], p1[1], 
                     v1[2], v1[1], 
                     vh1[2], vh1[1])

location$trustse <-  c(a1[2,2], NA, 
                       b1[2,2], NA,  
                       c1[2,2], NA,  
                       f1[2,2], NA,  
                       i1[2,2], NA, 
                       m1[2,2], NA,  
                       p1[2,2], NA,  
                       v1[2,2], NA,  
                       vh1[2,2], NA)

location$democracy <-  c(a2[2], a2[1],
                         b2[2], b2[1],
                         c2[2], c2[1],
                         f2[2], f2[1],
                         i2[2], i2[1],
                         m2[2], m2[1],
                         p2[2], p2[1],
                         v2[2], v2[1], 
                         vh2[2],vh2[1])

location$democracyse <-  c(a2[2,2], NA, 
                           b2[2,2], NA,  
                           c2[2,2], NA,  
                           f2[2,2], NA,  
                           i2[2,2], NA, 
                           m2[2,2], NA,  
                           p2[2,2], NA,  
                           v2[2,2], NA,  
                           vh2[2,2], NA)


location$platformmod <-  c(a3[2], a3[1],
                              b3[2], b3[1],
                              c3[2], c3[1],
                              f3[2], f3[1],
                              i3[2], i3[1],
                              m3[2], m3[1],
                              p3[2], p3[1],
                              v3[2], v3[1], 
                              vh3[2],vh3[1])

location$platformmodse <-  c(a3[2,2], NA, 
                                b3[2,2], NA,  
                                c3[2,2], NA,  
                                f3[2,2], NA,  
                                i3[2,2], NA, 
                                m3[2,2], NA,  
                                p3[2,2], NA,  
                                v3[2,2], NA,  
                                vh3[2,2], NA)

location$alliancemod <-  c(a4[2], a4[1],
                             b4[2], b4[1],
                             NA, NA,
                             f4[2], f4[1],
                             i4[2], i4[1],
                             m4[2], m4[1],
                             NA, NA,
                             v4[2], v4[1], 
                             vh4[2],vh4[1])

location$alliancemodse <-  c(a4[2,2], NA, 
                               b4[2,2], NA,  
                               NA, NA,  
                               f4[2,2], NA,  
                               i4[2,2], NA, 
                               m4[2,2], NA,  
                               NA, NA,  
                               v4[2,2], NA,  
                               vh4[2,2], NA)


location <- location %>% filter(beta == "Effect")

location <- merge(location,f, by = "Location")

#Plots
fa8a <- ggplot(na.omit(location), aes(x = reorder(Location, alliancemod),y= alliancemod #, shape = beta, color = beta
)) + 
  labs(#title = "Heterogeneity in campaign effects on alliance moderation",
       #subtitle = "By location (*Not measured in Planadas nor Coyaima)",
       x = "Location",
       y = "Effect size on behavioral moderation (ITT)") +
  geom_point()+   geom_pointrange(aes(ymin=alliancemod-alliancemodse, ymax=alliancemod+alliancemodse)) +
  scale_y_continuous(limits = c(-5, 8))+
  theme(text = element_text(size = 8))+
  theme(axis.text.x = element_text(face="bold", size=8,angle=40))+ geom_hline(yintercept=0, linetype="dotted", 
                                                                               color = "red", size=0.25)+
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) 


fa8b <- ggplot(location, aes(x = reorder(Location, platformmod),y= platformmod #, shape = beta, color = beta
)) + 
  labs(#title = "Heterogeneity in campaign effects on platform moderation",
       #subtitle = "By location",
       x = "Location",
       y = "Effect size on ideological moderation (ITT)") +
  geom_point()+   geom_pointrange(aes(ymin=platformmod-platformmodse, ymax=platformmod+platformmodse)) +
  scale_y_continuous(limits = c(-5, 8)) +
  theme(text = element_text(size = 8))+
  theme(axis.text.x = element_text(face="bold", size=8,angle=40))+ 
  geom_hline(yintercept=0, linetype="dotted", color = "red", size=0.25)+
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) 

fa8 <- ggarrange(fa8a, fa8b, legend = "none")

ggsave("fa8.tiff",plot = fa8, dpi = 300,  unit = "in")
ggsave("fa8.png",plot = fa8,unit = "in")


#### Figure A 9: Effect size and unit fractionalization by site #### 

coeff <- 1

fa9a <- ggplot(data = location, aes(x = reorder(Location, alliancemod)#, shape = beta, color = beta
)) + 
  geom_bar(aes(y= alliancemod), stat = "identity", fill = "lightgrey", color = "black") +
  geom_point(aes(y = uf*coeff), stat = "identity", color = "orange") + 
  scale_y_continuous(limits = c(-2, 5.5),
                     name = "Effect size on alliance moderation (ITT)",
                     sec.axis = sec_axis(~.*coeff, name="Fraciontalization")
  ) +
  labs(#title = "Effects on alliance moderation and site fractionalization",
    #subtitle = "By location",
    x = "Location")+
  theme(axis.title.y.right = element_text(color ="orange"))+
  theme(axis.text.x = element_text(face="bold", size=8,angle=40))+ 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) 


fa9b <- ggplot(data = location %>% filter(platformmod != "NA"), aes(x = reorder(Location, platformmod)#, shape = beta, color = beta
)) + 
  geom_bar(aes(y= platformmod), stat = "identity", fill = "lightgrey", color = "black") +
  geom_point(aes(y = uf/coeff), stat = "identity", color ="orange") + 
  scale_y_continuous(limits = c(-2, 5.5),
                     name = "Effect size on platform moderation (ITT)",
                     sec.axis = sec_axis(~.*coeff, name="Fraciontalization")
  ) +
  labs(#title = "Effects on platform moderation and site fractionalization*",
    #subtitle = "By location (*Not measured in Planadas nor Coyaima)",
    x = "Location")+
  theme(
    axis.title.y.right = element_text(color ="orange"))+
  theme(axis.text.x = element_text(face="bold", size=8,angle=40))+ 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) 

fa9 <- ggarrange(fa9a, fa9b, legend = "none")

ggsave("fa9.tiff",plot = fa9, dpi = 300,  unit = "in")
ggsave("fa9.png",plot = fa9,unit = "in")


#### Figure A 10: Effect size and workshop duration by sites #### 
load("rawdata_farcexperiment_english.RData")
farc$start_time <- trimws(farc$start_time)
farc$start_time <- gsub("a.m.", "AM",farc$start_time)
farc$start_time <- gsub("p.m.", "PM",farc$start_time)
farc$start_time <- gsub(": a. m.", " AM",farc$start_time)

farc$end_time <- trimws(farc$end_time)
farc$end_time <- gsub("a.m.", "AM",farc$end_time)
farc$end_time <- gsub("a.m.", "AM",farc$end_time)
farc$end_time <- gsub("p.m.", "PM",farc$end_time)
farc$end_time <- gsub(": a. m.", " AM",farc$end_time)

#Estimating duration of workshops
time <- dplyr::select(farc, municipio, start_time, end_time)
time$start_time <- (as.POSIXct(farc$start_time, format = "%H:%M"))
time$end_time <- (as.POSIXct(farc$end_time, format = "%H:%M"))
time$duration <- (time$end_time - time$start_time)
#Deleting workshops with incorrect time input (negative duration)
time$duration <- ifelse(time$duration < 0,NA,time$duration)
colnames(time)[1] <- "Location"

#Estimating average workshop duration and merging with location dataset

location <- merge(location, time %>% group_by(Location) %>% summarise(ws_duration = mean(duration, na.rm = T)), by = "Location")
coeff <- 60

fa10a <-  ggplot(data = location, aes(x = reorder(Location, alliancemod))) + 
  geom_bar(aes(y= alliancemod), stat = "identity", fill = "lightgrey", color = "black") +
  geom_point(aes(y = ws_duration/coeff), stat = "identity", color = rgb(0.2, 0.6, 0.9, 1)) + 
  scale_y_continuous(limits = c(-2, 5.5),
                     name = "Effect size on alliance moderation (ITT)",
                     sec.axis = sec_axis(~.*coeff, name="Duration (mins.)")) +
  labs(#title = "campaign effects on behavioral moderation and workshop duration",
       #subtitle = "By location",
       x = "Location")+
  theme(axis.title.y.right = element_text(color = rgb(0.2, 0.6, 0.9, 1)))+
  theme(axis.text.x = element_text(face="bold", size=8,angle=40))+ 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) 

fa10b <- ggplot(data = location, aes(x = reorder(Location, platformmod))) + 
  geom_bar(aes(y= alliancemod), stat = "identity", fill = "lightgrey", color = "black") +
  geom_point(aes(y = ws_duration/coeff), stat = "identity", color = rgb(0.2, 0.6, 0.9, 1)) + 
  scale_y_continuous(limits = c(-2, 5.5),
                     name = "Effect size on platform moderation (ITT)",
                     sec.axis = sec_axis(~.*coeff, name="Duration (mins.)")) +
  labs(#title = "campaign effects on behavioral moderation and workshop duration",
    #subtitle = "By location",
    x = "Location")+
  theme(axis.title.y.right = element_text(color = rgb(0.2, 0.6, 0.9, 1)))+
  theme(axis.text.x = element_text(face="bold", size=8,angle=40))+ 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) 


fa10 <- ggarrange(fa10a, fa10b, legend = "none")

ggsave("fa10.tiff",plot = fa10, dpi = 300,  unit = "in")
ggsave("fa10.png",plot = fa10,unit = "in")


#### Figure A 11: Distribution of education by site #### 
#Arauquita
arauquita <- ggplot(filter(imputed, municipio == "Arauquita"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) +
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Arauquita (N = 32)",
    x = "Schooling",
    y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
           color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Arauquita")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75) + 
  scale_x_continuous(breaks = 0:7)
  
ggsave("fa11a.png",plot = arauquita,unit = "in")

#Bogotá 
bogota <- ggplot(filter(imputed, municipio == "Bogotá"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey", binwidth = 1) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Bogotá (N = 16)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Bogotá")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11b.png",plot = bogota,unit = "in")


#Coyaima 
coyaima <- ggplot(filter(imputed, municipio == "Coyaima"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Coyaima (N = 29)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Coyaima")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7) + 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11c.png",plot = coyaima,unit = "in")


#Florencia
florencia <- ggplot(filter(imputed, municipio == "Florencia"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Florencia (N = 13)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Florencia")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11d.png",plot = florencia,unit = "in")


#Iconozo
iconozo <- ggplot(filter(imputed, municipio == "Iconozo"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Iconozo (N = 14)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Iconozo")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11e.png",plot = iconozo,unit = "in")


#Montañita 
montanita <- ggplot(filter(imputed, municipio == "Montañita"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Montañita (N = 32)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Montañita")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11f.png",plot = montanita,unit = "in")


#Planadas 
planadas <- ggplot(filter(imputed, municipio == "Planadas"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Planadas (N = 26)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Planadas")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11g.png",plot = planadas,unit = "in")


#Villavicencio 
villavicencio <- ggplot(filter(imputed, municipio == "Villavicencio"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Villavicencio (N = 22)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Villavicencio")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11h.png",plot = villavicencio,unit = "in")


#Vista Hermosa
vistahermosa <- ggplot(filter(imputed, municipio == "Vista Hermosa"), aes(x = schooling)) + 
  geom_histogram(position="identity", alpha = 0.5,colour = 1, fill = "grey",binwidth = 1 ) + 
  theme(panel.background = element_blank())+
  theme(axis.line = element_line(colour = "black")) +
  labs(title = "Vista Hermosa (N = 91)",
       x = "Schooling",
       y = "Frequency") +
  theme(text = element_text(size = 8))+
  geom_vline(xintercept = mean(imputed$schooling, na.rm = T),
             color = "purple", size=0.75)+
  geom_vline(xintercept = mean(filter(imputed, municipio == "Vista Hermosa")$schooling, na.rm = T), linetype = "dotted",
             color = "purple", size=0.75)+ 
  scale_x_continuous(breaks = 0:7, lim = c(0, 6)) 

ggsave("fa11i.png",plot = vistahermosa,unit = "in")



